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ADAPTIVE MULTIRESOLUTION SCHEMES WITH LOCAL TIME STEPPING FOR 
TWO-DIMENSIONAL DEGENERATE REACTION-DIFFUSION SYSTEMS 

MOSTAFA BENDAHMANE^, RAIMUND BURGER^, RICARDO RUIZ BAIER^, AND KAI SCHNEIDER^ 

Abstract. We present a fully adaptive multiresolution scheme for spatially two-dimensional, possibly de- 
generate reaction-diffusion systems, focusing on combustion models and models of pattern formation and 
chemotaxis in mathematical biology. Solutions of these equations in these applications exhibit steep gradi- 
ents, and in the degenerate case, sharp fronts and discontinuities. This calls for a concentration of compu- 
tational effort in zones of strong variation. 

The multiresolution scheme is based on finite volume discretizations with explicit time stepping. The 
' pJ multiresolution representation of the solution is stored in a graded tree ("quadtree"), whose leaves are the 

^^ non-uniform finite volumes on the borders of which the numerical divergence is evaluated. By a threshold- 

ing procedure, namely the elimination of leaves that are smaller than a threshold value, substantial data 
wA compression and CPU time reduction is attained. The threshold value is chosen optimally, in the sense that 

the total error of the adaptive scheme is of the same slope as that of the reference finite volume scheme. 

^^ Since chemical reactions involve a large range of temporal scales, but are spatially well localized (espe- 

"^^ cially in the combustion model), a locally varying adaptive time stepping strategy is applied. For scalar 

\^ ^ equations, this strategy has the advantage that consistence with a CFL condition is always enforced. Nu- 

• merical experiments with five different scenarios, in part with local time stepping, illustrate the effectiveness 

^I? of the adaptive multiresolution method. It turns out that local time stepping accelerates the adaptive 

^ multiresolution method by a factor of two, while the error remains controlled. 

>► 1. Introduction 

0\ 

l/^ Multiresolution techniques were first introduced by Harten [19] to improve the performance of schemes for 

CO one-dimensional conservation laws. Later on, these original ideas were extended to several kinds of related 

^"^ problems [ , , ], leading finally to the concept of fully adaptive multiresolution schemes [13, 14, 30, 35]. 

f"^ Overviews on multiresolution methods for conservation laws are given by Chiavassa, Donat, and Miiller [11] 

^D and Miiller [^^] . The basic aim of this approach is to accelerate a given finite volume scheme on a uniform 

25 g^id ^^ the cost of an at most controllable loss of accuracy, that is, the accelerated scheme should be of the 

same order than the original one. The principle of the multiresolution analysis is to represent a set of data 

given on a fine grid as values on a coarser grid plus a series of differences, called details^ at different levels 

of nested dyadic grids. These differences contain information on the local regularity of the solution. An 

^^ appealing feature of this data representation is that these details are small in regions where the solution is 

Cd smooth. By thresholding small details (cells whose coefficients are smaller than a prescribed tolerance are 

removed), a locally refined adaptive grid is defined. This threshold is chosen such that the discretization error 

of the reference scheme is balanced with the accumulated thresholding error introduced in each time step. 

Significant speed-up of the computation and data compression is achieved for long-time evolution problems, 

large systems, multidimensional domains, and solutions with sharp fronts. 
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The present paper serves two purposes. On one hand, the adaptive multiresolution scheme for par- 
abohc PDEs [ ] and strongly degenerate parabohc PDEs in one space dimension [ , ] is extended to 
two-dimensional systems of (possibly degenerate) parabolic PDEs. These equations produce solutions that 
vary smoothly wherever the solution causes the PDE to be parabolic, but produce sharp fronts, or even dis- 
continuities, close to solution values at which the equation degenerates, so adaptive multiresolution methods 
are a proper device to efficiently capture these fronts. Similar features (that is, solutions with steep gradi- 
ents) also appear in a combustion model of reaction-diffusion type. The analysis made in [ ] is extended 
here to the study of two interacting flame balls. In this case, an adaptive strategy is equally very useful, 
especially when the flame front is well localized in space, since fine grids are only needed in small subregions 
of the computational domain. We will utilize then an adaptive multiresolution scheme applied to a reference 
finite volume discretization with explicit time integration. 

On the other hand, chemical reactions are known to involve a large range of temporal scales, especially 
in long-time evolutions. Then an adaptive time stepping strategy is recommendable. Earlier efforts in this 
direction, which include [8, 12, 15, 18] and the references therein, were based on using the same time step 
to advance the solution on all parts of the computational domain, and controlling the time step through an 
embedded pair of Runge-Kutta schemes (known as Runge-Kutta-Fehlberg (RKF) schemes). In these proce- 
dures, one compares the numerical solution after each time step with an (approximate) reference solution, 
and adjusts the time step if the discrepancy is unacceptable. In contrast to this approach, we here adapt the 
locally varying time stepping strategy recently introduced for multiresolution schemes for conservation laws 
and multidimensional systems by Lamby, Miiller, and Stiriba [24] and Miiller and Stiriba [31]. This strategy 
is not precisely (time-) adaptive for scalar equations, since the time step for each level remains the same for 
all times. However, in the case of nonlinear systems, coupling of components entering the CFL condition 
makes it necessary to compute the time step after each iteration, according the evolving CFL condition, and 
therefore we have a scheme adaptive in time. Our results in terms of CPU time savings are encouraging 
and the strategy is consistent with a CFL condition, in contrast to the approach based one the RKF device. 
We mention that Miiller and Stiriba [ ] also combine local time stepping and multiresolution for implicit 
schemes, and that more details are also given in the germinal papers of Berger and Oliger [ ], Osher and 
Sanders [ ] and the references therein. 

The remainder of this paper is organized as follows. In Section 2 the reaction-diffusion systems studied 
herein are briefly described. In Section 3 we present the reference finite volume methods to which we apply 
the multiresolution device, and in Section 4 the adaptive multiresolution strategy, as well as the required 
graded tree data structure, are outlined. In Section 5 we analyze the error of the adaptive multiresolution 
scheme, and deduce the optimal choice of the threshold. This choice ensures that the discretization error of 
the reference scheme is balanced with the accumulated thresholding error which is introduced in each time 
step. In Section 6 we address the local time stepping strategy applied to the multiresolution strategy, and 
in Section 7, we outline the overall multiresolution procedure. Finally, in Section 8 the method is applied 
to different scenarios. Example 1 corresponds to a single-species reaction diffusion equation. Examples 2 
and 3 deal with the thermo-diffusive model for the interaction between flame balls. Examples 4 and 5 shows 
the results of Turing-type pattern formation produced by a reaction-diffusion system, and Example 6 arises 
from a model of chemotaxis with growth. Conclusions of our study are collected in Section 9. All numerical 
results clearly reveal high resolution and improvement in terms of compression of memory and savings in 
computational effort. 

2. A CLASS OF REACTION-DIFFUSION SYSTEMS 

2.1. A single-species reaction-diffusion model. Model 1 is the following initial-boundary value problem 
for a scalar reaction-diffusion equation, where x = (x^y) and {x^y^t) G Qt := Q x [0,T], ft C M?: 

Ut = f{u,^)^AA{u), (2.1a) 

7i(x,0) =iio(x) on(^, (2.1b) 

VA{u)'n = ouSt :=ai]x [0,T]. (2.1c) 
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This problem may serve as a scalar prototype degenerate reaction-diffusion model. Here, the zero-flux bound- 
ary condition (2.1c) implies that the reaction-diffusion domain is isolated from the external environment. 
For /(ix, x) = f{u), (2.1a) appears in [^^] in an ecological setting, where u denotes the population density 
of a species, and f{u) is its dynamics, where it is assumed that /(O) = and /'(O) 7^ 0. For example, 
f{u) = u{l — u) — V? /{I -^ V?) corresponds to the population dynamics of the spruce band- worm [ ], and 
models the growth of the population by a logistic expression and the rate of mortality due to predation by 
other animals. We modify this expression by a radial spatial factor, and use 

/(2i,x) := 10 ( exp(-5r)ii(l - ii) + (exp(-5r) - l) — ^ ) , r := a/(x - 0.5)^ + (^ - 0.5)^, (2.2) 

\ 1 + li^ y 

which means that the birth of individuals is concentrated near the center (0.5,0.5), and mortality increases 
with increasing distance from the origin. On the other hand, most standard spatial models of population 
dynamics simply assume that A{u) = Du^ where the constant diffusion coefficient D > measures the 
dispersal efficiency of the species under consideration. Motivated by Witelski [ ] , who advanced degenerate 
diffusion in the context of population dynamics, we utilize herein the strongly degenerate diffusion coefficient 

A{u) = l\ , ^°^"^"- (2.3) 

\D{u — Uc) otherwise, 

where Uc > is an assumed critical (threshold) value of u beyond which diffusion will take place. Model 1 
gives rise to Example 1 of Section 8. 

The difficulty in the well-posedness analysis of the problem (2.1) lies in the boundary condition (2.1c) 
when A is strongly degenerate. It is quite difficult to give a correct formulation of the zero flux boundary 
conditions. For the case of non-homogeneous Dirichlet boundary conditions, however, Mascia, Porretta, and 
Terracina [26] demonstrated existence and uniqueness of L^ entropy solutions. In the special case where 
the function A is strictly increasing, the classical framework of variational solutions of parabolic equations 
is sufficient to satisfy this wish. 

2.2. A two-species reaction-diffusion model. Model 2 is given by the following initial-boundary value 
problem for a reaction-diffusion system on Qt- 

ut = jf{u, v) + AA{u) on Qt, (2.4a) 

vt = jg{u, v) + dAB{v) on Qt, (2.4b) 

ii(x, 0) = iio(x), '^(x, 0) = 'Uo(x) for X G r^, (2.4c) 

VA{u) • n = \/B{u) • n = on St- (2.4d) 

We study this system in the context of two applications, namely as a model of combustion and as a two-species 
model of mathematical biology. 

For A{u) = B{u) = u^ d = 1/Le and 7 = 1, (2.4) represents a reduced dimensionless t her mo- diffusive 
model describing a combustion process, where Le is the Lewis number. We restrict ourselves to a simple 
chemical reaction with only two react ants and one product, the first reactant and the product being highly 
diluted in the second reactant; and we neglect gravity. Since the chemical reaction takes place in a lean 
premixed gas, we focus on the limiting reactant, and denote by v its normalized partial mass, while u 
represents normalized temperature. The reaction rates are given by an Arrhenius law: 

/,,,,:=^„«p( 'Szii^). ,(„,„):=-/(„,.), ,2.5, 

2 \<^(1 ~ ^) ~ ^ J 

where a and /3 are the temperature rate and the dimensionless activation energy, called Zeldovich number, 
respectively. In Example 2 of Section 8, this model is employed to simulate the interaction between two 
flame balls, as an extension of the applications of the same model that were considered in [ , ]. Here, a 
flame hall denotes a slowly propagating spherical flame structure in a premixed gaseous mixture. 
If radiation effects are taken into account, (2.4a) is replaced by 

Ut = -ff{u, v) + S{u) + AA{u) on Qt, (2.6) 
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where the dimensionless heat loss due to radiation S fohows the Stefan-Boltzmann law 

S{u) = p[{u^ a-^ - 1)^ - (a-i - 1)^] , (2.7) 

and the dimensionless coefficient p controls the radiation level. Conditions (2.4d) imply that the process 
takes place inside a box with adiabatic walls. See [ ] for details and a discussion of the case with one flame 
ball. The interaction of two flame balls including radiation is simulated in Example 3 of Section 8. 

On the other hand, (2.4) also arises in mathematical biology as a well-known reaction-diffusion system 
modelling the interaction between two chemical species with respective concentrations u and v. Under 
certain conditions, it produces stationary solutions with Turing-type spatial patterns [ 3, 40]. To simulate 
the formation of such a pattern, we here select the kinetics between each species due to Schnakenberg [38]: 

f{u^v) = a — u-\-u'^v^ g{u^v) = b — u^v. (2.8) 

Alternative choices of / and g that lead to Turing- type patterns are discussed in [32, 33]. For 

A{u) = B{u) = u, (2.9) 

this system has a uniform positive steady state {vP ^ v^) given by ix^ = a + 6 and v^ = b/{a-\-b)'^, where b > 
and a + 6 > 0, and under certain conditions, (2.4) has a unique solution. See for instance [ ] for the proof 
of existence and uniqueness. 

We recall from [33, Sect. 2.3] some results on the conditions under which (2.4) produces a diffusion-driven 
instability giving rise to Turing-type pattern in the non-degenerate case. A necessary condition is satisfaction 
of the inequalities fu^Qv < 0, fu9v-fvgu > 0, dfu^Qv > and {dfu^Qv)'^ -^d^fuQv- fvQu) > 0. Evaluating 
these inequalities for the system (2.4a), (2.4b) and the particular kinetics (2.8) yields the inequalities 

< 6 - a < (a + 6)^ (a + bf > 0, d{b -a)> {a^ 6)^ {d{b -a)-{a^ bff > 4d{a + bf. (2.10) 

To characterize the stationary pattern that arises from a choice of (a, 6) that satisfies (2.10), we define 

^ _ d{b -a)-{a^ b)^ ± ^/[d{b -a)-{a^ b)^]^ - 4d{a + b)^ 

The analysis of general rectangular domains [ ] implies that in the non-degenerate case, the unstable 
patterned solution of the initial-boundary value problem (2.4) is given by 

w(x, y, t) = 2^ ^nm exp(A(/c^)t) cos(n7rx) cos(m7r?/), 

where the constants c^m depend on a Fourier series of the initial conditions for w, and the summation takes 
place over all numbers n and m that satisfy 

jL~{a,b,d) =: kl < k^ = 7r^(n^ + m^) < kl \= -fL'^{a,b,d), 

and A(A:^) is the positive solution of the following equation, where fu^ fv^ 9u and gy are evaluated at (ix^, v^): 

X^ + X{e{l + d) - 7(A + gy)) + dk^ - ^{dfu + gv)k^ + l\fu9v - fv9u) = 0. 

Example 4 of Section 8 presents a numerical solution of (2.4) with the kinetics (2.8) and the diffusion 
coefficients (2.9), where parameters are chosen according to the preceding discussion such that indeed a 
Turing- type pattern is produced. On the other hand, in Example 5, we present a simulation where (2.9) is 
replaced by the degenerate diffusion functions 

A( \ /^ ioTu^u^, fo ioiu^v^, 

A(u) = i .1 . ' ^W = \ .1 . ^c,^c > 0. (2.11) 

\u — Uc otherwise \u — v^ otherwise. 

It turns out that even if the stability analysis does not apply to the degenerate case, our numerical experiments 
(Example 5) lead to the formation of a pattern. 

We mention that the mathematical analysis of the system (2.4) is still an open problem because of 
the presence of the boundary condition (2.4d). A successful technique for proving uniqueness of (entropy 
weak) solutions to degenerate parabolic equations with Dirichlet boundary condition is based on Kruzkov's 
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method [ ]. Related to this approach we mention that Holden, Karlsen, and Risebro [ ] prove existence 
and uniqueness of entropy solutions of weakly coupled systems of degenerate parabolic equations in an 
unbounded domain. Our system is an example of the degenerate reaction-diffusion system analyzed in [^^J, 
but is equipped here with the boundary condition (2.4d). 

2.3. A chemotaxis-growth system. We assume that 1] C M^ is convex, bounded and open. Model 3 is 
the following generalization of the Keller-Segel model [21, 22] for chemotactical movement: 

ut = \/ ■ {aVu - u\/x{v)) + g{u) on Qt, (2.12a) 

vt = h{u, v) + dAv on Qt, (2.12b) 

ii(x, 0) = 'Uo(x), 'u(x, 0) = 'Uo(x) on 1^, (2.12c) 

\/u-n = \/v-n = on St. (2.12d) 

The system (2.12) describes the aggregation of slime molds caused by their chemotactical features. Cell 
migration appears in numerous biological phenomena. In the case of chemotaxis, cells (or an organism) move 
in response to a chemical gradient. Specifically, (2.12) corresponds to the model proposed by Mimura and 
Tsujikawa [ ] for the spatio-temporal aggregation patterns shown by the bacteria Escherichia coli. This 
model incorporates four effects: diffusion, chemotaxis, production of chemical substance, and population 
growth. In the absence of growth, u = ii(x, t) usually represents the density of the cell population of the 
amoebae Dictyostelium discoideum, v = '^(x, t) is the concentration of the chemo- attract ant {cAMP: cyclic 
Adenosine Monophosphate)^ and x denotes the chemotactical sensitivity function, which may be given by 

X{v) = vv, v> 0, (2.13) 

where i^ is a chemotactical parameter. The function g takes into account the growth rate of the population, 
and can be given by 

g(u) = v?(\-u). (2.14) 

Moreover, a > and and d > are constant diffusion rates for both components. The function h describes 
the rates of production and degradation of the chemo- attractant; here, we choose 

h{u^v) = au — Pv^ a^P^O. (2.15) 

For this case it is known that if ^ i^o ^ L'^{Q), ^ vq G H^~^^{Q), and dQ is smooth enough, (2.12) 
possesses a unique global solution (see, e.g., [ ]); and if uq and vq are radial and ||'^o||li is sufficiently large, 
then ||i^(t) 11/^2 blows up in finite time. On the other hand, Efendiev, Klare, and Lasser [ ] analyzed the 
fractal dimension of the exponential attractor in dependence of u. Our Example 6 of Section 8 is based on 
examples presented in [ ], and presents numerical solutions of (2.12) for various values of u. 

3. Finite Volume Schemes 

We employ a standard finite volume scheme to discretize a reaction-diffusion equation, which is de- 
scribed here for a uniform grid. The rectangular spatial domain 1] C M^ is partitioned into control volumes 
(^ii)(i,j)GA, where A is an index set, defining a^:= [^^-1/2,^2+1/2] x [^^-1/2, ^j+1/2], Ax := x^+1/2 -^i-1/2, 
Ay := yj-\-i/2 — yj-i/2i foi" all {hj) ^ ^5 a^d Ax := min{Ax, Ay}. The cell average of a quantity q at time t 
is defined by 

3.1. Discretization of Models 1 and 2. The finite volume scheme is described here for (2.1) and as it 
applies to the first equation of (2.4); for the second equation of (2.4), we replace u by v., f{u^v) by g{u^v)^ 
and A{u) by dB{y). Integrating the respective equation and averaging over Vtij yields 

W\IL '"'^'^'^^'^''^WlIL ^("('^'*)'^^H'^'^)))^'^ +P^/{ /Mx,i))rfx, (3.2) 
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where V denotes the right-hand side of the PDE under consideration except for the reaction term. For the 
two-dimensional case and on a cartesian grid, V is discretized via 

The reaction term is approximated by fij ^ f{uij,Vij). If we incorporate a first-order Euler time integration 
for both components, then the corresponding interior marching formula for Model 2 is 

ul+' = ul + At7^- + MVij (5«,.), Ax), v^+^ = v^ + MrQij + ^At^i (5(^5), Ax), (3.3) 

where S{-) denotes the stencil utilized for computing Vij. According to [9, 20], this scheme is stable under 
the CFL condition 

A7(||/n||oo + ll/.lloo + ll^nlloo + \\9v\\oo) + 4/id(||A'|U + H^^loo) ^ 1- (3.4) 

Here A := At/ Ax, /i := At/Ax'^. 

3.2. Discretization of Model 3. We define the difference operators S^Vij := ±{Vi±ij — Vij) and S^Vij := 

=t(^i,j±i — Vij). Then a suitable second order difference operator for a general term V • {QVu) is 

V • {Q\/u) ^ ^^S^{Qi^i/2jS-Uij) + ^S:^{Qij^i/2S-Uij). 

Integrating the corresponding equations, averaging over ftij and discretizing yields the following interior 
marching formula: 

crAt.^__^ aAt^^^__^ At 



u 



n+l 



^r' =v^^^Ath{ul,v^^) + 0^+C^S- + ^,Sp;v-, 

1/2 := l{x\v^M^x\v^,^iK,^i). giVi/2,, := J (x'K)^i?, + x'(^IVl,.)^^l,,)• 



(3.5) 



Analogously to (3.4), the scheme (3.5) is stable under the corresponding CFL condition 

X{\\hu\\oo + ||/^.||oo + ll^'lloo) + 4/id(a + llx'lloo) ^ 1. (3.6) 

The left-hand sides of (3.4) and (3.6) obviously evolve in time, so in practice, at each time step we obtain 
At from these conditions, and A and /i are not constants; rather, they are adjusted in each time step. 

4. Conservative adaptive multiresolution discretization 

In this section we recall some basic properties of the multiresolution discretization and the data structure. 
For a more detailed description, we refer to [7, 35]. We will organize the numerical solution and corresponding 
differences at different levels, in a dynamic graded tree structure: whenever a node is included in the tree, 
all other nodes corresponding to the same spatial region in coarser resolutions are also included. The 
tree structure is mainly needed for ease of navigation, which contributes to accelerating the scheme; data 
compression would also be possible by other techniques. The adaptive grid corresponds to a set of nested 
dyadic grids generated by refining recursively a given cell depending on the local regularity of the solution. 

We denote by root the basis of the tree. In two space dimensions, a parent node has four sons, and the 
sons of the same parent are called brothers. A node without sons is called a leaf. A given node has s^ = 2 
nearest neighbors in each direction, called nearest cousins, needed for the computation of the fiuxes of leaves; 
if these nearest cousins do not exist, we create them as virtual leaves. Brothers are also considered nearest 
cousins. Figure 1 illustrates the graded tree structure. The leaves of the tree are the control volumes from 
which we form the adaptive mesh. Here, the property of the tree being graded means that grid refinement 
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Figure 1. Graded tree data structure (^'quadtree^'), after 



and coarsening is governed by that two neighboring control volumes cannot differ by more than one level 
in the tree. This is equivalent with the notion of "one-irregular rule" (see e.g., [ , ]). We denote by A 
the set of indices of existing nodes, by C{A) the restriction of A to the leaves, and by A^ the restriction of 
A to a multiresolution level /, < / < L. We denote by w^^^^^ = {uij^i^Vij^i) the vector of cell averages for 
both components of the solution (and the obvious simplification Wij^i = Uij^i for a single-species problem) 
located at spatial position (z, j) at level /, and by W^ the set of cell averages for all nodes at level /. To 
estimate the cell averages of a level / from those of the next finer level / + 1, we use the projection operator 
P^+i^/, which is exact, unique, and in our case is defined by 



w, 



ij^i 



1 
4 



7 . W2i+ei,2j+e2 

^i,e2G{0,l} 



,^+1- 



(4.1) 



To estimate the cell averages of a level / + 1 from those of level /, we employ the prediction operator P^^^+i, 
which provides an approximation w by interpolation of W/ at level / + 1. This operator is local in the sense 
that the interpolation for a son is made from the cell averages of its parent and the s nearest cousins of its 
parent; and it is consistent with the projection in the sense that it is conservative with respect to the coarse 
grid cell averages or equivalently, P^+i^z o P^^^+i = Id. For a regular grid structure in two dimensions, we 
use a polynomial interpolation introduced in [2]: 



where 



W2^+e„2,+e2,^+l = W,,,, , - (-l)^^Q. - {-l^Qy + {-l^'Q^y. ^i, 63 G {0, 1}, 



Qx •= 2^7n(Wi+n,j,/ — ^i-n,j,l), Qy •= /_^%{^i,j+p,l ~ ^i,j-p,l)^ 



n=l 



p=l 



Qxy •— 7 ^ 7n 7 ^lpV^i-\-nJ-\-p,l ~ ^^i-\-n,j-p,l ~ ^^i-nJ-\-p,l + ^ 



(4.2) 



(4.3) 



-P,i)' 



The chosen accuracy order of the multiresolution method for our cases isr = 5 + l = 3, where s is the 
number of required nearest uncles for each spatial direction. The corresponding coefficients are 7i = — ^ 
and 72 = jlg. Nevertheless, one may select here an arbitrarily higher order of accuracy. 

As stated before, the adaptive grid consists in the set of leaves >C(A), which forms a partition of the 
computational domain Q. 

A detail is the difference between the exact and the predicted value 



df 



'i,j,l '•~ ^ijyl ^i,3,l^ 



r1^ 



^^ij,l — '^ij.l- 



For multicomponent solutions, there are many possible definitions of a scalar detail dij^i that is calculated 
from the details of the components (in our case, d'^-^ 



and d'^-i), depending mainly on the nature of the 
problem. Roussel and Schneider [ ] utilize the Euclidian norm of the details, dij^i = {{df ■ ^)^ + (J^^- ^"i^^i^/^ 



""ij^l^ 



In our case, given the nature of our problems, we could simply select one component dij^i 

done by Sjogreen [ ] for the compressible Euler equations, but in order to guarantee that the computations 



df ■ 1, as was 
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Figure 2. Conservative flux computation for coarser levels. 

of the refinement and coarsening procedures are always on ttie safe side, in the sense that we always prefer 
to keep a node with a detail pair containing at least one value over the threshold 

ei = 22(^-^)s, (4.4) 

we will use dij^i = minjd^ • i-,d\ ■ i} and dij^i = max{(i^ • ^, (i^ • ^} for the refinement and coarsening procedures 
(see details in Algorithm 7.1), respectively, similar to Harten's choice [^^] for the Euler equations of gas 
dynamics. 

Since a parent has four sons, the consistency property of the prediction operator implies that knowledge 
of the cell average values of the four sons is equivalent to that of the cell average value of the parent node 
and three independent details. Repeating this operation recursively on L levels, we get the multiresolution 
transform on the cell average values M : W^ ^^ {Dl^ • • • , £^i, Wq), where Di = {dij^i)(^ij). This means that 
knowledge of the cell averages of all leaves is equivalent to that of the cell average value of the root and of 
the details of all nodes of the tree. 

After thresholding, i.e., removing nodes where the detail is below the prescribed tolerance \dij^i\ < si, 
where Si is given by (4.4); a safety zone is added to the tree, which means that one finer level is added to 
the tree in all possible nodes without violating the graded tree data structure. This is done by splitting each 
leaf into four new leaves in such a way that the new tree remains graded. This device, which was proposed 
e.g. in [1^, '^^J, ensures that the graded tree will represent adequately the solution in the next time step, 
and depends strongly on the assumption of finite propagation speed of sharp fronts. 

Now, assume that the tree has only two levels / and / + 1 (straightforwardly extensible to an arbitrarily 
larger tree). To ensure conservativity of the scheme, we compute only the fluxes at level / + 1 and we set 
the ingoing flux on the leaf at level / equal to the sum of the outgoing fluxes on the leaves of level / + 1 (see 
Figure 2) 

^(i+l,j,0^(i,i,0 = ^(2i+l,2j,^+l)^(2i+2,2j,^+l) + ^(2i+l,2i+l,^+l)^(2i+2,2j+l,^+l)- (4.5) 

This choice decreases the number of costly flux evaluations without loosing the conservativity in the flux 
computation, and this presents a real advantage when using a graded tree structure, see e.g. [ ]. This 
advantage is lost for a non-graded tree structure, in which case these data (fluxes for leaves on an immediately 
finer level) are not always available. 

It is important to remark that in the case of systems, since we are dealing with reaction-diffusion systems 
where species tend to attract each other, we manage the multiresolution framework and the data structure 
as one unified mesh with two components per control volume. This means that we construct only one 
graded tree and apply only one thresholding strategy for both species; however, there are other cases where 
it is preferable to organize the different species in separate adaptive meshes, for example when the species 
segregate spatially, as in systems of conservation laws modelling traffic flow and polydisperse sedimentation 



5. Error analysis of the adaptive multiresolution scheme 

Using the main properties of the reference finite volume scheme on a uniform grid at the finest level L, 
such as the contraction property in L^ norm, CFL stability condition and order of approximation in space, in 
[13, 35], the authors decompose the global error between the cell average values of the exact solution vector 
at the level L, denoted by w^^ = {u^^^v^^)^ and those of the multiresolution computation with a maximal 
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level L, denoted by w^j^, into two errors 

||w^x - WmrII < ||w^x - Wpvll + ||wfv - Wmr||. (5.1) 

The first error on the right-hand side, called discretization error ^ is the one of the reference finite volume 
scheme on a uniform grid at the finest level L. In many circumstances (see e.g. [^, 15]), this error can be 
bounded by 

II wi; - w^vll < C'i2-^^, Ci > 0, (5.2) 

where a is the convergence order of the finite volume scheme. Based on preliminary numerical experiments 
(obtained in a similar fashion as in [8]), for our examples we obtain the approximate value a = 2.18. 

For the second error, called perturbation error^ in [ ' ^] the authors assume that the details on a level / are 
deleted when they are smaller than a prescribed tolerance si. Under this assumption, they show that if the 
discrete time evolution operator is contractive in the chosen norm, and if the tolerance £i at the level / is 
given by (4.4), then the perturbation error accumulates in time and satisfies 

||wfv - w^rII < C2ne, C2 > 0, (5.3) 

where n denotes the number of time steps. At a fixed time T = nAt, this gives 

T 

||wfv - w^rII < ^2^e, C2 > 0. 

Motivated by their analysis, and according to the global CFL condition (3.4) of the reference finite volume 
scheme defined for the discretization of Model 2, we have 

At ^ 



_ ||/n||oo + ll/.lloo + IKIloo + Iklloo + AxM{\\A'\\^ + \\B'\\^) 

If we write Ax — ^J^\1~^ ^ this yields 



'll/«l|oo + ||M|oc + ||ff«||oc + ||ff.||oc + ^^2-i4d(||A'|U + ||B'||oo)' '" 

The main idea of the adaptive multiresolution scheme is to perturb the solution given by a finite volume 
scheme on a uniform discretization (reference mesh) in such a way that the total error, i.e., the error between 
the exact solution and the adaptive solution that is projected to the reference fine mesh, is of the same order 
as the discretization error. For this purpose, one has to balance the discretization error and the perturbation 
error. With this idea in mind, it is possible to derive in a similar fashion the optimal choice for the threshold 
parameter e for the adaptive multiresolution scheme. 

As in [^], we need that e oc 2~"^ or 

e|n|22i(||/„|U + IIMIoc + IKIIoc + Iklloo + v^2-^4d(||A'||oo + IIS'lloc)) ex 2-"^. 
In this way, for the numerical computations of Model 2, we may set the so-called reference tolerance sr to 

2-(a+2)L 
^^ ^^P.WuWoo + ll/.lloo + ll^nlloo + ll^.lloo) + W^^ ^d(\J^\U + Halloo)' ^^'^^ 

Analogously, for Model 3, the reference tolerance may be set to 

2-(ai+2)L 
^^ ^ ^ll^KII/lnlloo + ||/^.||oo + Iblloo) + Wm^M{a + llxlloo) ' ^^'^^ 

where olx is a value of the convergence rate for Model 3. 

Note that all the L^ norms in (5.4) and (5.5) are computed numerically. To determine an acceptable value 
for the factor C (which, of course, depends on T, Ci, C2 and C3), a series of computations with different 
tolerances are needed in each case, prior to final computations. Essentially, we select the largest available 
candidate value for C such that the same order of accuracy (same slopes for the error computation) as 
that of the reference finite volume scheme is maintained. This procedure basically generalizes the treatment 
in [^] of spatially one-dimensional strongly degenerate parabolic equations. In [ ] the authors prove for 
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scalar, one-dimensional, nonlinear conservation laws, that the threshold error is stable in the sense that the 
constant C is uniformly bounded and, in particular, does not depend on the threshold value ^r, the number 
of refinement levels L and the number of time steps n. In our case, even when a rigorous proof is still missing 
for the systems considered in the present work, from the previous deduction we see a similar behavior for C. 
We also mention that as in previous works [ , 8, 13], here the reference tolerance remains fixed for all 
times, though it is certainly possible to recompute the reference tolerance at each time step. 

6. Local time stepping 

We utilize a version of the locally varying time stepping strategy advanced by Miiller and Stiriba [31], and 
summarize here its principles. The basic idea is to enforce a local CFL condition by using the same CFL 
number for all levels, and the strategy consists in evolving all leaves on level / using the local time step 

Ati = 2^-^At, / = L-1,...,0, (6.1) 

where At = AtL corresponds to the time step on the finest level L. This strategy allows to increase the time 
step for the major part of the adaptive mesh without violating the CFL stability condition. 

Clearly, portions of the solution lying on different resolution levels need to be synchronized to obtain a 
conservative scheme. This will be achieved after 2^ time steps using Atf. all leaves forming the adaptive mesh 
are synchronized in time, so one time step with Ato is equivalent to 2^ intermediate time steps with AtL- 
In order to additionally save computational effort, we update the tree structure (refinement and coarsening) 
only each odd intermediate time step 1, 3, . . . , 2^ — 1 (as suggested in [l]), and furthermore, the projection 
and prediction operators are performed only on the levels occupied by the leaves of the current tree, i.e., we 
do not update the tree structure by prediction from the root cell, but from the coarsest leaves (we refer to 
this as partial grid adaptation). For the rest of the intermediate time steps, we use the current tree structure. 
Notice that the updating of the tree is still done in each global time step. For the sake of synchronization and 
conservativity of the fiux computation, for coarse levels (levels without leaves), we employ the same fiuxes 
{f>ij^i and fij^i) computed in the previous intermediate time step, because the cell averages on these levels 
are the same that in the previous intermediate time step. Only for finer levels (levels containing leaves), we 
compute fiuxes, and do so in the following way: if there is a leaf at the corresponding cell edge and at the 
same resolution level /, we simply perform a fiux computation using the brother leaves, and the virtual leaves 
at the same level if necessary; and if there is a leaf at the corresponding cell edge but on a finer resolution 
level / + 1 (in this case we refer to this edge as an interface edge)^ the fiux will be determined as in (4.5), that 
is, we compute the fiuxes at a level / + 1 on the same edge, and we set the ingoing fiux on the corresponding 
edge at level / equal to the sum of the outgoing fiuxes of the son cells of level / + 1 (for the same edge). We 
recall that the graded tree structure ensures that two neighboring control volumes of the adaptive mesh do 
not differ by more than one resolution level, which is equivalent to the satisfaction of the one-irregular rule. 
In order to always have at hand the computed fiuxes as in (4.5), we need to perform the locally varying 
time stepping recursively from fine to coarse levels. If at any instance of the procedure there is a missing 
value, we can project the value from the sons nodes or we can predict this value from the parent nodes. For 
illustrative purposes, we give an example of an interior first-order fiux calculation for Model 2, to complete 
a full macro time step, by the following algorithm (we show only the fiux calculation for (i, j, /) ^ (z + 1, j, /) 
and (i, j, /) -^ {i^j + 1, /) since the other fiuxes are obtained analogously): 

Algorithm 6.1 (Locally varying intermediate time stepping). 

(1) Grid adaptation (provided the former sets of leaves and virtual leaves). 

(2) do /c = 1, . . . , 2^ (and therefore the local time steps are n + 2~^, n + 2 • 2~^, n + 3 • 2~^, . . . , n + 1^ 

(a) Synchronization: 
do / = L, . . . , 1 

doi = l,...,|AJ,(0. j = l,...,|A|^(0 
if 1 ^ / ^ Ik-i then 

if (i^j^l) is a virtual leaf then 
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Zn+fc2-^ Zn+(/c-l)2-^ -n+/c2-^ n+(fc-l)2-^ 

Jij,l Jij,l ' ^ijV ^ij,/ 

endif 
else 

if (i^j^l) is a leaf then 

Zn+/c2-^ ^/ n+/c2-^ -n+fc2-^\ -n+/c2-^ ^^^-7^+^2-^ -n+fc2-^\ 

if (i + 1, j, /) is a leaf or {i^j + 1, /) is a leaf then compute fluxes by 

P(ij,i)^(ij+i,i) ' — Xm '^'^'^^''^'^^'^^ -^(^i,j,0) 

endif 

if (2i + 2, 2j, / + 1)^ (2z + 2, 2j + l, / + 1) are leaves (interface edges) then compute 
fluxes by 

Fii,j,l)^ii+l,j,l) ^- ^(2i+2,2j,^+l)^(2i+l,2j,^+l) + ^(2i+2,2j + l,^+l)^(2i+l,2j + l,^+l) 
F{ij^l)^{ij-\-l^l) ^- ^(2i,2j+2,Z+l)^(2i,2j + l,Z+l) + ^(2i+l,2j+2,Z+l)^(2i+l,2j + l,Z+l) 

endif 
endif 
endif 
enddo 
enddo 
(b) Time evolution: 

do / = 1, . . . , L, z = 1, . . . , \AUl), j = 1, . . . , |A|^(0 
if 1 ^ / ^ Ik-i then t/iere is no evolution: 

n+(/c+l)2-^ ,-.n+/c2-^ n+(fe+l)2-^ ,-.n+/c2-^ 

else 

Interior marching formula only for the leaves {i^j^ I): 

n+(fc+l)2~^ -n-\-k2~^ i ^.A-/- 7^n+/c2~^ i ^A-/-'n ( Q ( ^."^^+^2 



Vaa 7 ^ '^.•.- 



7M5i;r2" +dMA,45(f;,';+'=2" ),h{l)) 



endif 
enddo 

(c) Partial grid adaptation each odd intermediate time step: 

do / = L, . . . , 4 + 1 

Projection from the leaves. 
enddo 
do I = Iki ' ' ' iL 

Thresholding, prediction, and addition of the safety zone. 
enddo 
enddo 

Here, Ik denotes the coarsest level containing leaves in the intermediate step k (as introduced in [ ]), h{l) 
is the mesh size on level /, and |A|;2(/) is the size of the set formed by leaves and virtual leaves per resolution 
level / in the direction z. The interior marching formula is the modified marching formula (3.3) for Model 2, 
for the intermediate time steps /c = 1, . . . , 2^, for the leaf in the position (i, j) at level /. 

Our scheme is formulated for a first-order, explicit Euler time discretization. Generalizations for higher 
order schemes are given in [ ] for Crank-Nicolson schemes and in [ ] for Runge-Kutta schemes, respectively. 
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7. Algorithm implementation 
Now we give a brief description of the multiresolution procedure used to solve the test problems. 

Algorithm 7.1 (Multiresolution procedure). 

(1) Initialization of parameters. 

(2) Creation of the initial graded tree structure: 

(a) Create the root of the tree and compute its cell average value. 

(b) Split the cell, compute the cell average values in the sons and compute the corresponding details. 

(c) Apply the thresholding strategy for the splitting of the new sons: If dij^i > si then split the son 
(here we use di = mm{df^d]^}). 

(d) Repeat this until all sons have details below the required tolerance Si. 

(3) do n = 1, . . . , total dime steps 

(a) Determination of the leaves and virtual leaves sets. 

(b) Time evolution with global time step: Compute the discretized divergence operator for all the 
leaves. 

(c) Updating the tree structure: 

• Recalculate the values on the nodes and the virtual nodes by projection from the leaves. 
Compute the details for all positions (•, •, for I '^ h- 

if \dij^i\ < £i (here we use di = max((i^, (ij^)^ in a node and in its brothers then the cell 

and its brothers are deletable. 

endif 

• if some node and all its sons are deletable and the sons are leaves without virtual sons, 
then delete sons (coarsening). 

if this node has no sons and it is not deletable and it is not at level I = L, then 
create sons (refinement). 
endif 
endif 

• Update the values in the new sons by prediction operator from the former leaves. 
enddo 

(4) Output: Save meshes, leaves and cell averages. 

Here total dime steps stands for the total time steps needed to reach Tfinai using At as the maximum 
time step allowed by the CFL condition using the finest space step. 

When using a locally varying time stepping, replace step (3b) by the new step 

(3) do n = 1, . . . , total dim^e steps 

(a) Determination of the leaves and virtual leaves sets. 

(b) Time evolution with local time stepping: Compute the discretized divergence operator for all the 
leaves and virtual leaves 

(c) do /c = 1, . . . , 2^ (k counts intermediate time steps) 

• Compute the intermediate time steps depending on the position of the corresponding leaf 
as explained in Section 6. 

• if /c is odd then update the tree structure: 

— Recalculate the values on the nodes and the virtual nodes by projection from the 
leaves. Compute the details in the whole tree. 

if |^i,j,z| < Si (here we use di = max((i^, d^)^ in a node and in its brothers then the 

cell and its brothers are deletable. 

endif 

— if some node and all its sons are deletable and the sons are leaves without virtual 
sons, then delete sons (coarsening). 

if this node has no sons and it is not deletable and it is not at level I = L, then 
create sons (refinement). 
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endif 
endif 

— Update the values in the new sons by prediction operator from the former leaves. 
endif 
enddo 

(Now, after 2^ intermediate steps, all nodes are synchronized.) 
enddo 

Here total dime _steps for the new step stands for the total time steps needed to reach Tfinab with Ato as the 
maximum time step ahowed by the CFL condition using the coarsest space step. 

Notice that with such a process, we obtain high-order approximation in the smooth regions and mesh 
refinement near discontinuities as a consequence of the polynomial accuracy in the multiresolution prediction 
operator, even if the reference finite volume scheme is low-order accurate. 

Bihari and Harten [^] use the following quantity, which we call data compression rate, 

„.= ^^^^^y (71) 

^" 7VL.iVi,/22i + |£(A)|' "•'■ ' 

to measure the possible improvement in data compression, whose feasibility in turn strongly depends on a 
smart implementation to navigate inside the tree, see for example [ ]. Here, N^^^Ly is the number of cells 
in the finest grid, and |>C(A)| is the size of the set of leaves. We measure the speed-up between the CPU time 
of the numerical solution obtained by the FV method and the CPU time of the numerical solution obtained 
by the multiresolution method: V = (CPUtime)Fv/(CPUtime)MR- 

To measure errors between a reference solution u and an approximate solution umr obtained using 
multiresolution, we will use L^-errors: Cp = \\u^ — i^mrIIp' P = I7 2, 00, where 

600 = max r^7 7 ~ ^mr? i tA^ 

Here ^^^Ri 7 l ^^ ^^^ value on the finest level L obtained by prediction from the corresponding leaf. 

8. Numerical results 

8.1. Example 1: Single-species model. For this example, we consider (2.1) with a strongly degenerate 
diffusion term (2.3), where we choose D \= \ and u^ := 0.5, a square domain ^ = [0, 1]^, and the function 
f{u,x.) given by (2.2). Figure 3 displays the numerical solution starting from the smoothly varying initial 
function 

uo{x, y) = 0.5(1 + sin(l.l(x — cos{0.7y))) cos(0.5{y — sin(1.3x))). 

We choose a maximal resolution level of Nl = 256^ = 65536 control volumes on the finest grid. Figure 4 
illustrates how the factor C in (5.4) is selected in our case as the optimal value from a finite selection of test 
values (each value giving a different value for the reference tolerance (5.4) sr). Figures 4 (a) and (b) indicate 
that for all displayed levels, the multiresolution procedure is in every case (for different values of C) cheaper 
(in terms of both acceleration and memory savings) than the corresponding reference FV computation on 
the finest grid. Figure 4 (c) indicates that the computations obtained using C = 1.0 x 10^ (and hence 
5r = 9.43 X 10~^) are sufficiently accurate, in the sense that with these choices, we keep the same slope for 
the I/-*^ -error as the FV calculations while increasing V and 77. We remark that here, as in previous works 
that use similar methods (see e.g. [ ]), there actually exists a range of threshold parameters that preserve 
the same slope for errors with respect to reference solution, for which C = 1.0 x 10^ is an average value. 
Here, we compute errors using a reference FV solution on a fine grid with Nl = 2048^ = 4194304 control 
volumes. (Here and in all other examples, we calculate errors in the approximate sense with respect to a 
reference solution.) 
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Figure 3. Example 1 (single- species model): Numerical solution (left) and leaves (right) of 
the corresponding tree data structure at times t = (top), t = 0.5 (middle) and t = 3 
(bottom). 



8.2. Examples 2 and 3: Interaction between two flame balls. We study (2.4) as a dimensionless model 
for the interaction between two flame balls of different sizes. We consider a square domain Q = [—30,30]^ 
and that the walls are sufficiently far from the flame balls so that their influence is negligible. Physical 
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Figure 4. Example 1 (single- species model): (a) data compression rate rj, (h) speed-up 
factor V and (c) L^ -errors for different levels L and values of C. The simulated time is 
t = 0. 



parameters characterizing the gaseous mixture and the combustion process are chosen as in [36, 37]. We use 
the parameters a = 0.64 and /3 = 10. The initial data is given by u{x^ y^ 0) = 1^0(^17^2), v{x^ y^ 0) = '^0(^1, ^2) 



with rf = {x — xi) 



•r 



{x — X2) + y , where 



^0(^1,^2) := 



max{exp(l — ri/a),exp(l — r2/6)} 



if ri < a or r2 < 6, 
otherwise, 



^0(^,^2) := 1 -^0(^1,^2). (8.1) 



In Example 2, we simulate the process without radiation, i.e., p = and hence S{u) = 0. We set the Lewis 
number to Le = 1. Here xi = —7.5, X2 = 7.5 and a = 1.8, 6 = 2.5 are the respective x— position and initial 
radii of the two flame balls. This choice ensures that there is no interaction between the two flame balls at 
t = and that there is no extinction of the flame balls. We simulate the process until t = 10, and Figure 5 
shows from left to right the temperature and reaction rate configuration obtained using the fully adaptive 
multiresolution scheme, and the position of the dynamic graded tree leaves, which form the corresponding 
adaptive mesh. The different times correspond from the top to the bottom to: before (t = 2), during {t = 4) 
and after (t = 10) direct interaction between the two flame balls, when the balls tend to create a new 
circular flame structure. We choose the following multiresolution parameters: the maximal resolution level 
L = 10 corresponding to Nl = 512^ = 262144 control volumes in the finest grid, and the reference tolerance 
Sr = 4.94 X 10-^. 

For comparison purposes, we introduce the global chemical reaction rate 



R{t)'= [[ f{u,v)dxdy. 
J Jn 



Errors in different norms, reaction rates, information on data compression and speed-up rate for different 
methods at three different times are depicted in Table 1. Due to the particular shape of solutions, which is 
nearly constant away from the combustion front, by using multiresolution, one can obtain very high rates of 
data compression, speed-up and low errors. 

In Example 3, we simulate the case with radiation, i.e. we use (2.6) and (2.7), where p = 0.05 and 
Le = 0.3. Now xi = —5, X2 = 5 and a = 0.5, 6=1 are the respective x— position and initial radii of 
the two flame balls. We simulate the process until t = 10 and Figures 6 and 7 show the scenario for this 
case. First, the balls grow spherically, tend to create a new flame structure, and then their fronts tend to 
extinguish when they touch each other, while the radiation effect causes the entire flame front to split. Here 
the maximal resolution level is set to L = 10 corresponding to Nl = 512^ = 262144 control volumes in the 
finest grid, and the reference tolerance is Er = 7.43 x 10~^. 
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MR 
FV 
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98.7942 



Table 1. Example 2 (interaction of two flame halls without radiation): Corresponding sim- 
ulated time, speed-up rate V , compression rate rj, errors, and total reaction rate R{t) for the 
u species. 



Notice that the multiresolution procedure automaticahy detects the higher gradient regions and uses this 
information to adaptively represent the solution by the refinement and coarsening of the mesh, i.e., by the 
adaptive addition and removal of control volumes on these areas. 

The L^, L^ and L'^ errors between the numerical solution obtained by our multiresolution scheme for 
different multiresolution levels L and the reference solution (obtained by finite volume approximation in a 
uniform fine grid with 2^'^^ control volumes) for Example 2 are depicted in Figure 8 (c) and (d). The slopes 
indicate a rate of convergence slightly better than two. 

For Example 3, we apply the locally varying time stepping (LTS) strategy detailed in Section 6. We choose 
the maximum CFL number allowed by (3.4), which is CFLq = 1 for the coarsest level. For the remaining 
levels we use CFL^ = 2^ CFLq, which means that we perform each macro time step with At = Ato = 2^AtL 
as given by (6.1). 

In Figure 9 we compare speed-up, data compression rate and total reaction rate for the finite volume 
reference scheme, the multiresolution scheme with global time step, and the multiresolution method with 
level-dependent time stepping. Notice that with LTS, the speed-up rate is approximately doubled for all 
times, while the compression rate and the total reaction rate remain of the same order as the multiresolution 
computation with global time step. 

8.3. Examples 4 and 5: A Turing model of pattern formation. We select the parameters a = 
—0.5, b = 1.9, d = 4.8 and 7 = 210. According to the discussion of Section 2.2, these parameters allow 
diffusion-driven instabilities to evolve. The initial concentration distribution is a normally distributed random 
perturbation around the stationary state (u^^v^) for the non-degenerate case, with a variance lower than 
the amplitude of the final patterns, see Figure 10. For the case of non-degenerate diffusion (Example 4), we 
use A{u) and B{u) as given by (2.9). For these parameters, the steady state is {u^ = 1.4,1'^ = 0.96939). 

In Example 4 we choose a maximal resolution level of Nl = 256^ = 65536 control volumes in the finest 
grid and a reference tolerance given by sr = 2.6 x 10~^. The time step is the maximum allowed by the 
CFL condition (3.4). Table 2 summarizes the speed-up rate, compression rate and errors in different norms 
between the numerical solution by multiresolution and the fine-mesh finite volume reference solution for 
different times. We depict errors between our multiresolution scheme and a reference FV solution with 
Nl = 1024^ = 1048576 control volumes in the finest grid, for different multiresolution levels L in Figure 12 
(c) and (d). In this case, the slopes equally indicate a rate of convergence slightly larger than two. Concerning 
the computation of errors, in Examples 4, 5 and 6 the system is evolved until the "random noise", which 
is imposed as an initial condition on the finest grid, has been smoothed sufficiently; then, this solution is 
projected on coarser levels to obtain auxiliary initial conditions for all the needed levels. 

For Example 5, we use the degenerate diffusion coefficients (2.11) with Uc = 1.2 and Vc = 0.7, and employ 
again the kinetics (2.8), but this time we choose the parameters a = —0.5, b = 1.9, d = 4.8 and 7 = 395. 
We select a maximal resolution level of Nl = 256^ = 65536 control volumes in the finest grid, with a 
reference tolerance given by s^i = 3.59 x 10~^. From Table 3 we see that the multiresolution algorithm 
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Figure 5. Example 2 (interaction of two flame balls without radiation): Numerical solution 
for species u (left) and leaves of the corresponding tree data structure (right) at times t = 2 
(top), t = A (middle) and t = 10 (bottom). 



allows significant acceleration and data compression rate are significantly increased by the multiresolution 
algorithm with very good accuracy. Figure 13 indicates that due to the degeneracy of the diffusion given by 



18 



BENDAHMANE, BURGER, RUIZ, AND SCHNEIDER 



u(x,y,t) 




10.9 
0.8 



Position of the leaves 



0.7 
-0.6 
0.5 
0.4 
0.3 
0.2 
0.1 




u(x,y,t) 




■ 0.9 
0.8 



pQSitJQn of ftie le^rvBE 



-0.7 
-0.6 
0.5 
0.4 



0.3 

10.2 
0.1 




Figure 6. Example 3 (interaction of two flame halls with radiation): Numerical solution 
for species u (left) and and leaves of the corresponding tree (right) at times t = (top) and 
t = 5 (bottom). 



Time 


V 


V 


Species 


L-^— error 


L^— error 


L"^— error 


t = 0.05 


7.16 


11.3783 


u 

V 


6.81 X 10-"^ 
4.09 X 10-"^ 


4.76 X 10-^ 
3.92 X 10-"^ 


3.46 X 10-3 
5.38 X 10-"^ 


t = 0.25 


9.29 


11.9756 


u 

V 


8.37 X 10-"^ 
4.22 X 10-^ 


6.94 X 10-^ 
5.43 X 10-"^ 


9.93 X 10-3 
8.48 X 10-"^ 


t = 1.50 


11.87 


14.4739 


u 

V 


9.26 X 10-"^ 
4.30 X 10-"^ 


2.71 X 10-"^ 
9.77 X 10-^ 


2.44 X 10-2 
8.39 X 10-3 



Table 2. Example 4 (Model 2 with non- degenerate diffusion): Corresponding simulated 
time, CPU ratio V , compression rate r] and componentwise errors. 



(2.11), and in contrast to Example 4, species u exhibits patterns with steeper gradients, and especiahy at 
t = 0.25 and t = 1.5, singularities appear. 
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Figure 7. Example 3 (interaction of two flame halls with radiation): Numerical solution 
for species u and leaves of the corresponding tree (right) at times t = 12 (top) and t = 20 
(bottom). 



Time 


V 


V 


Species 


L-^— error 


L^— error 


L"^— error 


t = 0.10 


6.32 


12.5438 


u 

V 


6.31 X 10-"^ 
4.98 X 10-"^ 


5.82 X 10-"^ 
5.37 X 10-"^ 


2.72 X 10-3 
9.46 X 10-"^ 


t = 0.25 


9.79 


10.3457 


u 

V 


6.12 X 10-"^ 
3.91 X 10-^ 


2.46 X 10-^ 
9.22 X 10-"^ 


3.03 X 10-3 
9.92 X 10-"^ 


t = 1.50 


11.60 


10.1984 


u 

V 


3.42 X 10-"^ 
2.63 X 10-"^ 


7.34 X 10-"^ 
4.98 X 10-"^ 


3.40 X 10-3 
2.81 X 10-3 



Table 3. Example 5 (Model 2 with degenerate diffusion): Corresponding simulated time, 
-up rate V , compression rate rj and componentwise errors. 



8.4. Example 6: Chemotaxis-growth system. For Example 6, in (2.12) we consider a square domain 
Q = [0, 16]^ and fix the parameters a = 0.0625 and d = 1. The function h{u, v) is given by (2.15) with a = 1 
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Figure 8. Example 2 (interaction of two flame balls without radiation): (a) speed-up rate V , 
(h) data compression rate rj, for different levels, at time t = 4.0; (c) errors \\umr — '^refill; 

II^MR - ^reflb; ||^MR - ^ref||oo and (d) \\vmK " ^ref||l; ||^MR - ^reflb and \\vmK " ^ref ||oo 

respectively for different levels L, at time t = A. 



and /3 = 32. The growth function g{u) for the species u is given by (2.14), and the chemotactical sensitivity 
is given by (2.13). This configuration corresponds to the model of chemotaxis and growth presented in 
[27], which is further analyzed in [ ]. Similarly to [ ], the initial datum is {uq^vq) = (1 + £(x),l/32), 
where ^(x) is a particular smooth perturbation which goes to zero near (8,8). We simulate the process 
until the solution reaches inhomogeneous stationary states, and we present three cases corresponding to 
different values of v^ which is responsible for the complexity of the spatial patterns. For example, for i^ = 7 
Figure 14 (middle) shows labyrinth-shaped patterns and for z/ = 10 (bottom), single filaments and spots. 
The corresponding adaptive meshes were generated with Nl = 512^ = 262144 control volumes in the finest 
grid, with £r = 8.43 x 10~^. For all these cases we implement locally varying time stepping, so we will choose 
the maximum CFL number allowed by (3.6), CFLq = 1 for the coarsest level and CFL^ = 2^CFLo for finer 
levels. From Figure 15 we can observe that if we incorporate the local time stepping strategy, a substantial 
gain (a factor slightly lower than 2, which is consistent with the results by Lamby, Miiller, and Stiriba [ 4]) 
is obtained in speed-up rate when comparing with a multiresolution calculation using global time stepping. 
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Figure 9. Example 3 (interaction of two flame balls with radiation): Time evolution of 
speed-up, data compression, and total reaction rates; and L^ -errors for different methods. 
L = 10 multiresolution levels and reference tolerance Sr = 7.43 x 10~^. 



The errors are computed from a reference solution in a grid with Nl = 2048^ = 4194304 control volumes. 
We conclude that the errors are kept of the same slope that the errors obtained with a global time step. 

The compression rate rj for both methods is lower than in the previous examples, which could be explained 
by the complexity and density of the spatial patterns in this particular example. 



9. Conclusions 

This paper describes an adaptive multiresolution scheme combined with a locally varying time stepping 
used to approximate solutions of a class of two-dimensional reaction-diffusion systems in Cartesian geometry. 
Several numerical examples show that the adaptive multiresolution mechanism with a suitable choice of the 
threshold value represents a gain in CPU time while the errors are kept of the same order as the reference 
finite volume method. In Examples 3 and 6, we also see that the local time stepping strategy is responsible 
for a gain in CPU time speed-up for a factor of about 2. Also, the errors between the solution using local 
time stepping and a reference solution are of the same order that the solution obtained by the adaptive 
multiresolution with global time stepping. 
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Figure 10. Examples 4 (^nd 5 (Turing pattern-formation): Initial data uo{x^y) (left) and 
vo{x,y) (right). 



The motivation to employ explicit schemes only is the following. Even though implicit methods allow 
larger time steps, we need to iteratively solve a nonlinear system in each time step, using e.g. Newton- 
Raphson method. The number of iterations is usually controlled by measuring the residual error, and 
cannot be controlled a priori. Thus, it appears difficult to assess the true benefits of a time-stepping strategy 
if the basic time discretization is an implicit one. 

On the other hand, of course, for the Turing-type pattern formation problem, it is conceded that patterns 
appear when one eigenvalue goes from negative to positive. At steady state (when the pattern is visible) all 
the eigenvalues again have negative real part. Thus to converge to steady state once the domain of attraction 
of the pattern is reached, implicit methods offer significant advantages since they can use larger and larger 
time steps. 

We remark that for hyperbolic problems, the incorporation of an implicit time discretization to the MR- 
LTS strategy can possibly form a substantial improvement in the speed-up rate, as presented in [31]. 
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